Design of a multi-epitope vaccine candidate against Brucella melitensis

Brucella is a typical facultative intracellular bacterium that can cause zoonotic infections. For Brucella, it is difficult to eliminate with current medical treatment. Therefore, a multi-epitope vaccine (MEV) should be designed to prevent Brucella infection. For this purpose, we applied the reverse vaccinology approach from Omp10, Omp25, Omp31 and BtpB. Finally, we obtained 13 cytotoxic T lymphocyte (CTL) epitopes, 17 helper T lymphocyte (HTL) epitopes, 9 linear B cell epitopes, and 2 conformational B cell epitopes for further study. To keep the protein folded normally, we linked AAY, GPGPG, and KK to CTL epitopes, HTL epitopes, and B cell epitopes, respectively. The N-terminal of the vaccine peptide is supplemented with appropriate adjuvants to enhance immunogenicity. To evaluate its immunogenicity, stability, safety, and feasibility, a final MEV containing 806 amino acids was constructed by linking linkers and adjuvants. In addition, molecular docking and molecular dynamics simulations were performed to verify the affinity and stability of the MEV-TLR4. Then, codon adaptation and in silico cloning studies were carried out to identify the possible codons for expressing the MEV. In animal experiments, the results demonstrated that the MEV had high immunogenicity. Collectively, this study provided a theoretical basis for the development of a Brucella vaccine.

www.nature.com/scientificreports/ Prediction of T-cell epitopes of proteins. T-cell epitopes are provided by class I (MHC I) and II (MHC II) MHC molecules that are recognized by two distinct subsets of T-cells, CD8, and CD4 T-cells, respectively 19,20 . Depending on the polymorphism of the HLA allele, the T cell response is different. We selected the alleles with high frequency in Xinjiang (HLA-A*1101, HLA-A*0201, HLA-A*0301, HLA-DRB1*0701, HLA-DRB1*1501, and HLA-DRB1*0301) to predict the cytotoxic T lymphocyte (CTL) epitopes and helper T lymphocyte (HTL) epitopes 21 . IEDB (http:// tools. immun eepit ope. org/) and NetCTLpan1.1 server (https:// servi ces. healt htech. dtu. dk/ servi ce. php? NetCT Lpan-1.1) were applied to predict the CTL epitopes of the target protein 22,23 . During CTL epitope prediction 1 was set as a threshold value for epitope identification. Since NetCTLpan1.1 counts from 0, we add one to the starting sequence number of all NetCTLpan1.1 to facilitate comparison with the results of IEDB. Accordingly, we listed the top 10 high-score epitopes for each prediction software. IEDB and NetMHC-IIpan-4.0 (https:// servi ces. healt htech. dtu. dk/ servi ce. php? NetMH CIIpan-4.0) were applied to predict the HTL epitopes of the target protein 24 . When predicting HTL epitopes, the default threshold for NetMHC-IIpan4.0 is set to the original threshold. Finally, we selected two software overlapping sequences as T cell dominant epitopes for the proteins.
Prediction of B-cell epitopes of proteins. B-cell epitopes are the antigen portion binding to the immunoglobulin or antibody. B-cell epitopes include linear epitopes and conformational epitopes. BCPREDS (http:// ailab. ist. psu. edu/ bcpred/ predi ct. html) was used to predict the selection of dominant B-cell linear epitopes with a specificity of 75%. Then, discontinuous (conformational) epitopes were performed by Ellipro (http:// tools. iedb. org/ ellip ro/) of the IEDB. After prediction, we also deleted the sequences with less than 5 amino acids and selected the sequences with high scores and appropriate lengths. In the end, linear epitopes and conformational epitopes are considered for the design of MEV.
The interaction between HLA alleles and T-cell epitopes: molecular docking. As common interactions of the HLA alleles, HLA class I (HLA-A*02:01) and HLA class II (HLA-DRB1*01:01) were selected to look for structural associations between the HLA alleles and the T cell epitopes. The molecular docking was performed in the HDOCK server.

Development of a novel MEV.
Non-allergic and non-toxic epitopes were selected for the MEV. Then, the selected epitopes were joined together with appropriate linkers to construct the vaccine sequence 25 . Different linkers were used to link B-cell epitopes together, the KK linker was added; for linking HTL epitopes together, the GPGPG linker was appended; and concurrently, for connecting CTL epitopes, the AAY linker was used. Human β-defensin-3 sequence (Serial No. Q5U7J2) and PADRE sequences were linked with the help of EAAAK linkers at the N-terminus. Finally, a polyhistidine tag was added at the C terminus to obtain the complete vaccine protein sequence 26 .
Evaluation of vaccine construct. Bioinformatics software was used to synthesize and analyze the properties of the MEV. First of all, ProtParam software was applied to analyze the physicochemical properties of the MEV. VaxiJen was applied to analyze the antigenicity with an accuracy of 70 to 89%. Then, SOLpro (http:// scrat ch. prote omics. ics. uci. edu/) was used to predict the solubility of MEV. The server has an overall accuracy of 74.15% and the threshold is 0.5 27 . Finally, AllergenFP was used to analyze its allergenicity.
Prediction of secondary and tertiary structure. To predict secondary structures of MEV, the SOPMA (http:// npsa-pbil. ibcp. fr/ cgibin/ npsa_ autom at. pl? page=/ NPSA/ npsa_ sopma. html) secondary structure prediction tool was used. RoseTTAFold uses a machine learning-based method and includes a three-track network to process sequence, the distance between atoms, and coordinate information to predict the tertiary structure of the MEV 28 .
Assessment of model quality. The quality of the tertiary structure models of the MEV was evaluated using the SWISS-MODEL Structure Assessment (https:// swiss model. expasy. org/ assess) and ProSA-web server (https:// prosa. servi ces. came. sbg. ac. at/ prosa. php) 29 . The assessment of the model quality is based on the Z-score and the dotted structure in the Ramachandran plot.

Molecular docking.
When the MEV enters the organism, the body recognizes it as a synthetic antigen and activates the immune system 30 . Therefore, we performed molecular docking to show the binding between proteins and immune molecules. We use the HDOCK server (http:// hdock. phys. hust. edu. cn/) to perform the molecular docking between the MEV and the TLR4 (PDB ID: 4G8A) immune receptor 31 . Key interacting residues were deduced from the protein-ligand interaction diagram generated by LigPlot+ in combination and their 3D structure was visualized by PyMOL.
Immune simulation. The vaccine enters the body as an antigen and also causes a corresponding immune response. C-ImmSim (https:// kraken. iac. rm. cnr. it/C-IMMSIM/) was applied to simulate the extent and type of immune responses induced by the MEV in humans 32 . In the population of the Xinjiang region, high-frequency alleles of HLA-A*1101, HLA-A*0201, HLA-B*5101, HLA-B*3501, HLA-DRB1*0701, and HLA-DRB1*1501 were selected for analysis 21 . Simultaneously, the server simulates three sections that represent three different anatomical regions in mammals, including the bone marrow, the thymus, and the lymph. And three different intervals were 1, 84 and 168 33  www.nature.com/scientificreports/ tion parameter random seed was set to 12,345, the simulation volume was set to be 50, and the simulation steps were set to 1050.

Molecular dynamics simulation. Molecular dynamics (MD) simulation studies of the MEV-TLR4 com-
plex were performed in Gromacs 2021.3 34 . Being solvated in a rectangular box of TIP3P waters, the resulting systems were extended up to a minimum cutoff of 15 Å from the protein boundary. Being added to the protein surface, Cl− or Na+ ions were used to neutralize the total charges of the systems. Afterward, the Amber ff14SB force field was employed for the protein in all of the MD simulations 35 . Using combined steepest descent and conjugate gradient method, the initial structures were fully minimized. Under canonical ensemble for 0.2 ns, the systems were then gently annealed from 10 to 300 K with a weak restraint of 15 kcal/mol/Å. What's more, 100 ps of density equilibration was performed under an isothermal-isobaric ensemble at a target temperature of 300 K and the target pressure of 1.0 atm using Langevin-thermostat and Berendsen barostat with a collision frequency of 0.002 ns and pressure-relaxation time of 0.001 ns. Finally, a productive MD run of 10,000 ps was performed for all the complex systems after proper minimizations and equilibrations.
Optimization of codons and in silico cloning. The XHOI and BamHI restriction endonuclease sites were selected and the codons of the vaccine were analyzed and optimized using the JAVA codon adaptation tool (http:// www. jcat. de/). The DNA sequence of the plasmid was obtained from the Plasmid Files, and pET-28a (+) was selected as the vector. Then SnapGene 5.3.2 was used to analyze the suitable endonuclease sites in the multiple cloning site (MCS) region, and the target gene sequence of the MEV was inserted into the plasmid to complete the silicon cloning. And the MEV was synthesized and purified by SynPeptide Co Ltd (Shanghai, China). Lastly, the MEV was detected as endotoxin-negative by Limulus amebocyte lysate (LAL) before use.
Simulation polymerase chain reaction and agarose gel electrophoresis. We  The mice were randomly divided into two groups (N = 10 per group), the control group, and the MEV group. Mice were instilled intranasally (i.n, 15 μL per nostril) with PBS (Control group), and 30 μg multi-epitope vaccine (MEV group) respectively for 15 days. The mice were executed and the respiratory lymph nodes (head, neck, and lower respiratory lymph node) and spleen were removed under sterile conditions. Then, the homogenate was prepared for cell isolation. The trachea was cannulated, and the lungs were gently lavaged 3 times with 1 mL of sterile PBS to obtain bronchoalveolar lavage fluid. Using a Pasteur capillary pipette without anticoagulant, blood was collected from the retro-orbital sinus and placed in coagulation microtubes. Finally, the blood was centrifuged and the serum was collected.
ELISA. To evaluate the immunogenicity of MEV, antigen-specific IgG1 and IgG2a titers in the immunized mice were tested by enzyme-linked immunosorbent assay (ELISA www.nature.com/scientificreports/ respiratory lymph nodes cells and splenocytes were added in triplicate to wells at 1 × 10 5 cells/well in a final volume of 100 μL. Fifthly, cells were stimulated with the heat-killed RB51 (10 μg/mL) (a source of Brucella protein without LPS), and synthetic Toxoplasma SAG4 peptide (10 μg/mL). ConA (10 μg/mL) at 50 μL/well and PBS at 50 μL/well were used as the positive and negative controls, respectively. After that, the plates were incubated for 24 h at 37 °C in an atmosphere containing 5% CO2. After washing, 50 μL of biotin-labeled secondary antibody (BD PharMigen, San Diego, CA) was added and the plates were incubated at room temperature for 2 h. And the plates were washed with PBST and 50 μL streptavidin alkaline phosphatase (Amersham Life Science, Australia) was added and incubated at room temperature for 2 h. By adding BCIP/NBT developer (Moss, Inc, Pasadena, MD, USA) to each well for incubation for 40 min at room temperature, the spots were developed. Finally, the plates were washed in water to stop the reaction and Spot forming units (SFU) in each well were counted using an ELISPOT Bio Reader-4000 (BIOSYS, GmbH, Germany). The results were expressed as IFN-γ spot forming cells (SFC) per million cells.
Approval for animal experiments. The animal study was reviewed and approved by The Animal Experiment Medical Ethics Committee of the First Affiliated Hospital of Xinjiang Medical University. We confirmed that all experiments were performed in accordance with relevant named guidelines and regulations. We confirmed that the authors complied with the ARRIVE guidelines. Sequence retrieval. The MAFFT in Jalview with the highest accuracy ( Fig. 2) showed that these four proteins have several homologies. Thereinto, the high homology of Omp25 and Omp31 indicates that these proteins originated from the same gene during the evolution of B. melitensis, and suggests that they may play a similar role in the immune response.
Prediction of T-cell epitopes of proteins. First of all, epitopes with high scores for each software were selected. Secondly, VaxiJen was applied to analyze the antigenicity of the epitopes. Furthermore, AllergenFP v1.0 (http:// www. ddg-pharm fac. net/ Aller genFP/) was used to detect whether epitopes are allergens. Then, the ToxinPred (https:// webs. iiitd. edu. in/ ragha va/ toxin pred/ design. php) was applied to predict the toxicity of the epitopes 37 . Finally, 13 CTL dominant epitopes and 17 HTL dominant epitopes were obtained ( Table 1). The antigenicity of these epitopes has higher antigenicity, and they were non-toxicity and non-allergenicity.
Prediction of B-cell epitopes of proteins. Through bioinformatics analysis, we obtained 9 dominant linear B-cell epitopes. In the analysis of B-cell conformational epitopes, we obtained 5 conformational epitopes of Omp10, 4 conformational epitopes of Omp25, 8 conformational epitopes of Omp31, and 5 conformational epitopes of BtpB. Finally, we selected non-toxicity and non-allergenicity B cell conformational epitopes (Fig. 4) to construct the vaccine ( Table 2).
The interaction between HLA alleles and T-cell epitopes: molecular docking. To evaluate the structural associations with HLA allele(s), we performed molecular docking simulations to find the interaction between the HLA alleles and the T cell epitopes. The interaction between HLA-I and CTL epitopes showed that the docking score was − 232.24 and ligand RMSD was 47.94 Å. The interaction between HLA-II and HTL epitopes showed that the docking score was − 326.30 and ligand RMSD was 143.19 Å. All the results indicated that docked complexes had a good affinity (Fig. 5).
Development of a novel multi-epitope vaccine. The MEV construct contained the dominant epitopes.
There were 13 CTL epitopes, 17 Th epitopes, 9 LB epitopes, and 2 CB epitopes in the vaccine. All the dominant epitopes for the construction of MEV were shown in Fig. 6. The number of vaccine residues of the MEV was 806 amino acids.  As illustrated in the figure, "A" stands for Omp2a, "B" stands for Omp2b, "C" stands for Omp10, "D" stands for Omp16, "E" stands for Omp19, "F" stands for Omp22, "G" stands for Omp25, "H" stands for Omp28, "I" stands for Omp31, "J" stands for BtpB, "K" stands for BLS, "L" stands for BCSP31, and "M" stands for P39.         44.17% random coil, and 26.55% extended strand. The proportion of alphahelix, β-turn, random coil, and the extended strand was consistent with the tertiary structure. This indirectly indicated that the prediction of the tertiary structure was reasonable. Then, we showed the PDB format of the MEV in Discover Studio, where pink in the model is more likely to represent the donor and green to represent the receptor (Fig. 7).
Assessment of model quality. ProSA-web was applied to assess the quality of the model, and the Z-score in the software indicated the quality of the model (Fig. 8). The light blue region indicated the structural group with the source of X-rays, and the dark blue region indicated the structural group with the source of nuclear magnetic resonance (NMR). The Z-score of MEV was − 7.32. The Z-score within the range characteristic for native proteins is indicative of a correct structure. To ensure the accuracy of the structure, we further evaluated the quality of the model using SWISS-MODEL's structure evaluation service, which is called the Ramachandran plots. The dark green region in the Ramachandran plots indicated the allowed region, the light green region in the diagram indicated the maximum allowed region, and the blank region in the diagram indicated the disallowed region. It can be seen that most of the amino acids (i.e., dotted structures) fall almost completely inside the interval range. In conclusion, the tertiary structure of the MEV is the correct model with high confidence.

Molecular docking. Molecular docking could effectively predict whether the ligand and the receptor could
interact through the principle of energy minimization and the complementarity of the spatial structure in the region of the receptor active site. We obtained many clusters by molecular docking, of which we selected the number one cluster for presentation. The appropriate molecular docking results showed that the docking score was -328.00 and the ligand RMSD was 145.30 Å. After selecting the best docking structure (Fig. 9A), 3D structures were then visualized using PyMOL (Fig. 9B) and a 2D interaction map was visualized using LigPlot + (Fig. 9C) 38 .
Immune simulation. The C-ImmSim server was employed to predict the stimulated immune response against the MEV through a computational approach 32 . The simulator works by Miyazawa and Jernigan proteinprotein potential measurements, applying input amino acid sequences to represent pathogen and lymphocyte receptors, and assessing the level of immune response elicited when the MEV is injected into the body. B cells are mainly involved in stimulating the immune response and cause an increase after each MEV injection. The number of B cells peaked in vivo after the last injection (Fig. 10A). Th cells and TC cells are important components of the T cell subpopulation. After three doses of the vaccine, the growth trend of Th cells (Fig. 10B) was less than that of B cells but reached a peak after three doses. In contrast, TC cells (Fig. 10C), NK cells (Fig. 10D), DC cells (Fig. 10E) and EP (Fig. 10F) generally remained relatively stable with dynamic changes. IgM and IgG concentrations are continually elevated in the primary and secondary immune responses accompanying antigen reduction (Fig. 10G). Finally, MEV can also elicit a highly responsive cytokine response (Fig. 10H), causing significant elevations in IFN-γ, TGF-β, IL-10 and IL-18.

Molecular dynamics simulation.
MD simulations of the MEV with TLR4 structure need to be conducted. For this, a 10 ns simulation of the MEV-TLR4 was performed, using the GROMACS and MDrun tools. The structural parameters pertaining to the MD simulations like the root mean square deviations (RMSD), the root mean square fluctuations (RMSF), the radius of gyration (RoG), and the potential energy analysis of the MEV-TLR4 were estimated to assess their conformational and structural stability (Fig. 11). The structural variations of the MEV-TLR4 were examined by RMSD values from 0 to 10 ns, showing a steady increase at the begin-   The compactness of the complex structure was determined by enumerating the Rog data of the complexes from the X, Y, and Z-axis. And the fluctuation range was less than 0.05 nm, which indicated that the backbone atoms in the complex structure are relatively stable in terms of gyration. Non-bonded interaction energy was calculated that are two types of short-range potential: Lennard-Jones short-range (LJ-SR) and Coulombic short-range (Coul-SR) potential. We calculated the total sum of LJ-SR and Coul-SR potential of ligands. For 10 ns simulation, the total residue interaction energy result was − 734.058 kJ/ mol. The obtained sum of LJ-SR and Coul-SR potential were in line with the experimental activity of each ligand.

Optimization of codons and in silico cloning. The codon adaptation index (CAI) is shown for cod-
ing sequences designed to maximize optimal codon usage (CAI = 1 indicates that the sequence was designed to have a codon adaptation index as close as possible to 1.00). The optimized CAI value of the MEV is 0.96. The ideal percentage range of GC content of the gene was between 30 and 70%, after optimization, the average GC content was changed to 70%. The BamHI and XHOI restriction enzyme sites were inserted at the N and C ends of the optimized codons, respectively. Then, the target gene was inserted into the MCS domain. Finally, specific restriction enzyme cutting sites were considered (Fig. 12).
Simulation polymerase chain reaction and agarose gel electrophoresis. According to the above principles, we designed the forward primer (5′-CTT ATC TCG AGG AGG CCG -3′) with a length of 18, Tm value of 56.31, and a GC content of 60%; the reverse primer (5′-CAT GGA TCC GTG GTGGT-3′) with a length of 17, Tm value of 56.05 and the GC content of 58%. BLAST results showed that the MEV shares the highest homology with the Homo sapiens zona pellucida glycoprotein 2 (ZP2). Afterward, the target gene of the MEV was amplified using SnapGene. We carried out the simulated agarose electrophoresis at a concentration of 1% and selected TBE with the better buffering ability of the target gene, vector, and recombinant plasmid in buffer solution. Finally, the amount of DNA was consistent with previous predictions. The sequence of the target gene was 2418 bp, the sequence of pET-28a (+) was 5639 bp, and the sequence of recombinant plasmid was 7753 bp (Fig. 13).

The results of ELISA and ELISPOT experiments.
The MEV has to be tested in an animal model to see if the humoral and cellular immune response can be elicited. ELISA was used to verify whether humoral immunity could be produced. The results showed that the MEV group produced significantly higher levels of IgG1 represents the X-axis direction, "blue" represents the Y-axis direction, and "green" represents the Z-axis direction. (D) Potential energy analysis that included Lennard-Jones short-range and Coulombic short-range. www.nature.com/scientificreports/ and IgG2a (Fig. 14A) than the control group, which was dominated by the elevation of IgG2a (Fig. 14B). And the difference was statistically significant (P < 0.01). In the mucosal immune system, the results showed that the MEV group produced higher levels of sIgA (Fig. 14C) than the control group, and the difference was statistically significant (P < 0.01). In the same way, ELISPOT was used to verify whether cellular immunity could be produced. The results (Fig. 14D,F) showed that the MEV group produced significantly higher levels of IFN-γ than the control group, and the difference was statistically significant (P < 0.01). There were no statistically significant differences between the PBS group and the SAG4 peptide group (P > 0.05). The number of IFN-γ-producing cells (visualized as spots) can show the difference between results more visually (Fig. 14E,G). Our study demonstrated that MEV can elicit good humoral and cellular immunity in an animal model.

Discussion
Brucellosis is a reemerged zoonotic infectious disease with complex clinical symptoms and high mortality rates 39 . The vaccine is still preferred for brucellosis prevention, but there are no existed Brucella vaccines for humans 40,41 . Several latest types of research have been reported recently, which encouraged the reverse vaccinology top down approach to get dominant epitopes from the whole proteome of viral and bacterial [42][43][44] . In our study, MEV was innovatively designed and its immunogenicity was evaluated in animal models. Proteins with high antigenicity are possible targets for the construction of MEVs 45 . Therefore, based on the current research, we analyzed 13 major Brucella-associated proteins. After that, we chose Omp10, Omp25, Omp31, and BtpB with high antigenicity, stability, and hydrophilicity to build the MEV. At present, based on these four candidate proteins, there are no relevant reports for the construction of MEVs. In the current study, Omp25 of B. melitensis was shown to have significant homology with Omp31 of B. melitensis. At length, we again  www.nature.com/scientificreports/ confirmed the high homology of these proteins when comparing the sequences, which provides a theoretical basis for future mechanisms of brucellosis pathogenesis 46 . The signal peptide influenced the start of protein translation, and the different primary structures of the signal peptides even affected the folding and transport of the proteins 47 . Therefore, we removed the signal peptide sequences of Omp10 (1-31), Omp25 , and Omp31 (1-31), respectively. However, BtpB does not possess the signal peptide, so we considered the whole amino acid sequence can be further analyzed.
The objective of predicting HTL epitopes and CTL epitopes is to find the short peptide sequence within an antigen that stimulates CD4+ or CD8+ T cells in vivo 48 . Furthermore, the prediction of B cell epitopes of pathogenic organisms broadens the types of antigenic peptides and helps to improve the immunogenicity of antigens 49 . Therefore, four proteins were analyzed using IEDB and NetCTLpan1.1 server, which obtained 13 dominant CTL epitopes. Using IEDB and NetMHCIIpan4.0, 17 dominant HTL epitopes were acquired from four proteins. Using BCPREDS, 9 dominant linear B-cell epitopes were obtained from four proteins. Afterward, 2 dominant conformational B-cell epitopes from four proteins were obtained using IEDB. Ultimately, docked complexes between HLA alleles and T-cell epitopes had a good affinity.
It has shown that single-epitope-based vaccines had weak immunogenicity in the present study, but an MEV that consists of multiple epitopes can induce the activation of both cellular and humoral immune responses and enhance its immunogenicity 50 . Therefore, an effective MEV should be designed to include effective epitopes that can induce HTL, CTL, and B cells 51 . Three types of linkers that can keep proteins in are intrinsically ordered and normally folded states, including flexible linkers, rigid linkers, and cleavable (or self-cleavable) linkers 52,53 . The vaccine sequence has been constructed by combining CTL epitopes, HTL epitopes, and B-cell epitopes with the The number of T cells that can produce IFN-γ in splenocytes specimens after antigen stimulation. (G) The representative ELISPOT spot diagram. One-way ANOVA was used for the Statistical Analysis. Data are shown as mean ± SD. P < 0.05 was considered statistically significant (**P < 0.01). www.nature.com/scientificreports/ AAY, GPGPG, and KK linkers, respectively. After that, we constructed a vaccine with multiple epitopes. Adjuvants maintain the stability of the peptides and increase their immunogenicity 54 . There are several adjuvants used in reverse vaccinology studies, such as human β-defensin-3 (hBD3) and PADRE sequence 55,56 . HBD3 is a cationic peptide with immunomodulatory effects on both innate and acquired immune responses 55 . The PADRE sequence can reduce the effect of human HLA-DR polymorphism, and enhance the long-term immune response by inducing CD4+ T cells 57 . Finally, a complete MEV was obtained by adding adjuvants to the peptide sequence and adding the histidine sequence at the end. In the early stages of vaccine design, some features of the MEV should be considered. In reverse vaccinology, theoretically, the molecular weight of the protein used to design the vaccine should be less than 110 KD 58 . The results of this project revealed that the MEV with 806 amino acids in length and molecular weight of 86 KD was a soluble protein. In addition, it had an antigenicity of 0.9750 (greater than the threshold) and had no allergenicity. These results indicated that the MEV had good stability, hydrophilicity, antigenicity, solubility, and non-allergenicity. After that, the prediction of the secondary structure showed that the β-turn and random coil accounted for 8.44% and 44.17%, respectively. The high proportion of β-turn and random coil in the MEV suggests that the protein is likely to form antigenic epitopes 59 . The RoseTTAFold had ultra-high tertiary structure prediction accuracy by successively integrating and converting the three-track network information of one-dimensional sequence, two-dimensional distance, and three-dimensional coordinates. Then, we used RoseTTAFold to predict the tertiary structures of four proteins and the MEV. Finally, the ProSA-web and the structure assessment service of SWISS-MODEL were used to verify the tertiary structure quality of the MEV. The results showed that the tertiary structure of MEV has high accuracy and a high approximation coefficient with the correct structure. In general, the tertiary structure of the MEV also exhibited a high number of β-turn and random coil, which is consistent with the results predicted by the secondary structure, indicating that the MEV has efficient antigen potential.
TLR4 is a transmembrane receptor protein with extracellular leucine-rich repeated domains and a cytoplasmic signaling domain, specifically recognizing endogenous molecules released from damaged or ischemic tissues termed pathogen-associated molecular patterns (PAMPs), and damage-associated molecular patterns (DAMPs) 60 . TLR4 recognizes lipopolysaccharides (LPS) of gram-negative bacteria and acts as a receptor for molecular docking with the MEV. Through docking, it can be seen that there is a strong interaction force between molecules through the atomic interaction interface between TLR4-MEV. The results suggest that the vaccine can develop stable connections with immune receptors so that it can be transported throughout the host body 61 .
To understand and explore the immune response profile and immunogenicity, an immune simulation study was conducted 62 . By carrying out the immune simulation experiment, we found that B and T cells in the body gradually increased after three doses of the vaccine and peaked at the time of the last dose of the vaccine. In addition, the MEV could also increase the levels of IFN-γ, TGF-β, IL-10 and IL-18. In particular, IFN-γ plays a major role in cell-mediated immunity by enhancing the bactericidal responses of macrophages, promoting the differentiation of Th1 cells, inducing B cell antibody class switching, enhancing cytotoxic responses of NK cells, and stimulating antigen presentation to T cells 63 . High levels of immune cells and cytokines indicate that the vaccine has great potential to induce a good immune response in humans. To further analyze the structural stability of the docked TLR4-MEV, a 10 ns MD simulation of the complex was performed. Taking into account the starting structure as reference for the TLR4-MEV, RMSD, RMSF, and ROG were computed to evaluate TLR4-MEV conformation changes and their stability. All the results suggested that the TLR4-MEV has strong stability.
The expression of the MEV in a suitable E. coli expression vector is pivotal for the production of recombinant proteins 64 . JCat was used for the optimization of codons, to increase the expression of MEV in the E. coli. CAI of optimized vaccine protein was 0.96 and the average GC content of our sequence was 70% showing good expression possibility of vaccine protein in E. coli host. After amplification, we obtained a 7753 bp recombinant plasmid that provided a favorable theoretical basis for further research.
To evaluate the effects of the MEV, we chose the specific antibodies as reliable indicators for humoral immunity and used IFN-γ as representative indices of cellular immunity. Higher IgG1, IgG2a, and sIgA demonstrated that MEV can elicit good humoral immunity in animal models. Simultaneously, higher levels of IFN-γ was demonstrated that MEV can elicit good cellular immunity in animal model.
Structure-based reverse vaccinology is an effective approach using rational vaccine design, but the resulting immune response does not represent effectiveness against all strains 65 . Because of this, we chose B. melitensis for further study, which is often infected by humans. In the present study, an MEV including CTL epitopes, HTL epitopes, and B cell epitopes, coupled with appropriate adjuvants, was designed to enhance the immune response in animal models. The results showed that this MEV becomes an excellent and suitable candidate against B. melitensis.

Conclusion
Brucella is a gram-negative coccobacillus that can cause zoonotic diseases. Reverse vaccinology presented in the present work may produce new knowledge about vaccines against B. melitensis. Dominant epitopes were fused using linkers and adjuvant to enhance the MEV's immunogenicity. Physiochemical properties, antigenicity, allergenicity, as well as solubility and tertiary structure analysis, of MEV were found to be very satisfactory. Moreover, the docked complexes during the simulation revealed that a strong and stable binding interaction of MEV with TLR4. In addition, this study used the C-ImmSim server to define the MEV, which is highly immunogenic. Finally, the study demonstrated that MEV can elicit both humoral and cellular immune responses in animal models.

Data availability
The data was derived from public domain information: Uniprot database (https:// www. unipr ot. org/) and PDB library (https:// www. rcsb. org/). The data that support the findings of this study are available in the methods. The data that support the findings of this study are available from the corresponding author upon request. There are no restrictions on data availability. If you have any questions, please contact corresponding author. www.nature.com/scientificreports/ Publisher's note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.